knitr::opts_chunk$set(echo = TRUE)
Sys.setenv('MAPBOX_TOKEN' = 'pk.eyJ1IjoiYWxvbnNvMjUiLCJhIjoiY2tveGJseXJpMGNmcDJ3cDhicmZwYmY3MiJ9.SbThU_R8YGE1Zll-nNrZKA')

library(readxl) # lectura archivos excel.
library(tidyverse);library(corrplot) # librearias de tratamiento de data
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6      ✔ purrr   0.3.4 
## ✔ tibble  3.1.8      ✔ dplyr   1.0.10
## ✔ tidyr   1.2.0      ✔ stringr 1.4.1 
## ✔ readr   2.1.2      ✔ forcats 0.5.2 
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## corrplot 0.92 loaded
library(sf); library(raster); library(rgdal) # librerias espaciales basicas
## Linking to GEOS 3.9.1, GDAL 3.4.3, PROJ 7.2.1; sf_use_s2() is TRUE
## Loading required package: sp
## 
## Attaching package: 'raster'
## 
## The following object is masked from 'package:dplyr':
## 
##     select
## 
## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## 
## rgdal: version: 1.5-32, (SVN revision 1176)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.4.3, released 2022/04/22
## Path to GDAL shared files: C:/Users/Alonso/AppData/Local/R/win-library/4.2/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/Alonso/AppData/Local/R/win-library/4.2/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.5-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
library(mapview); library(mapedit); library(mapdeck) # librerias espaciales visualizacion especifica
## Warning: package 'mapdeck' was built under R version 4.2.2
## 
## Attaching package: 'mapdeck'
## 
## The following object is masked from 'package:tibble':
## 
##     add_column
library(leaflet); library(leaflet.extras)
library(htmltools)
#ms <- mapdeck_style("satellite")
mapviewOptions(basemaps = "Esri.WorldImagery")

Publicación Internacional:

Certificado de Aceptacion:

Study to Determine Levels of Cadmiumin Cocoa Crops Applied to Inland Areas of Peru

Study to Determine Levels of Cadmiumin Cocoa Crops Applied to Inland Areas of Peru

Paper:

Acreditation

Acreditation

Video de Presentación

El link es el siguiente: ECITEC UNI 2020

Análisis de la Información Recolectada:

Lectura de Data

Cacao <- read_xlsx(path = "data_IIC1/BD_Cd.xlsx", col_names = TRUE)
str(Cacao)
## tibble [24 × 7] (S3: tbl_df/tbl/data.frame)
##  $ Xeast         : num [1:24] 518642 518644 518690 518701 518706 ...
##  $ Ynorth        : num [1:24] 9032291 9032343 9032334 9032232 9032284 ...
##  $ Cdsoil_mgkg   : num [1:24] 1.75 1.6 1.9 1.7 1.25 1.75 1.5 1.4 1.4 1.35 ...
##  $ Pbsoil_mgkg   : num [1:24] 5.41 5.03 7.24 5.81 5.91 5.62 6.57 5.26 5.01 7.02 ...
##  $ Znsoil_mgkg   : num [1:24] 70.5 99.4 77.1 97.1 74.7 ...
##  $ pHapprox      : num [1:24] 7.06 6.55 6.28 8 6.7 6.4 6.53 5.79 5.41 6.39 ...
##  $ Ceapprox_uS_cm: num [1:24] 39.9 10.3 5 126.7 8 ...
head(Cacao, n=5)
tail(Cacao, n=5)
summary(Cacao)
##      Xeast            Ynorth         Cdsoil_mgkg     Pbsoil_mgkg   
##  Min.   :518642   Min.   :9031868   Min.   :1.100   Min.   :4.740  
##  1st Qu.:519033   1st Qu.:9031913   1st Qu.:1.350   1st Qu.:5.247  
##  Median :519183   Median :9032036   Median :1.575   Median :5.850  
##  Mean   :519108   Mean   :9032052   Mean   :1.634   Mean   :5.929  
##  3rd Qu.:519221   3rd Qu.:9032134   3rd Qu.:1.863   3rd Qu.:6.525  
##  Max.   :519556   Max.   :9032343   Max.   :2.550   Max.   :7.240  
##   Znsoil_mgkg       pHapprox     Ceapprox_uS_cm  
##  Min.   :11.80   Min.   :5.410   Min.   :  4.50  
##  1st Qu.:44.79   1st Qu.:6.160   1st Qu.:  9.00  
##  Median :58.70   Median :6.375   Median : 13.50  
##  Mean   :58.66   Mean   :6.494   Mean   : 26.11  
##  3rd Qu.:71.89   3rd Qu.:6.588   3rd Qu.: 40.02  
##  Max.   :99.40   Max.   :8.000   Max.   :126.69
# Verificar nulos
Vacios <- sapply(Cacao, function(x) sum(is.na(x)))
Vacios
##          Xeast         Ynorth    Cdsoil_mgkg    Pbsoil_mgkg    Znsoil_mgkg 
##              0              0              0              0              0 
##       pHapprox Ceapprox_uS_cm 
##              0              0
#Con mapa:
library(Amelia)
## Loading required package: Rcpp
## ## 
## ## Amelia II: Multiple Imputation
## ## (Version 1.8.0, built: 2021-05-26)
## ## Copyright (C) 2005-2022 James Honaker, Gary King and Matthew Blackwell
## ## Refer to http://gking.harvard.edu/amelia/ for more information
## ##
missmap(Cacao,col=c("black","grey"),legend = FALSE)
## Warning: Unknown or uninitialised column: `arguments`.
## Unknown or uninitialised column: `arguments`.
## Warning: Unknown or uninitialised column: `imputations`.

# Verificar alguna condicion y eliminar valores
# id <- which(is.na(Cacao$pHapprox)) # colocar cualquier condicion
# Cacao <- Cacao[-id, ]

id <- which(Cacao$Cdsoil_mgkg>= 1.4 & Cacao$Cdsoil_mgkg <= mean(Cacao$Cdsoil_mgkg, na.rm=TRUE))
Cacaonew <- Cacao[id, ] 
Cacaonew %>% dplyr::select(Pbsoil_mgkg) %>% summarise(media_Pb = mean(Pbsoil_mgkg))

Exploracion Basica

Univariante

par(mfrow=c(2,2))
hist(Cacao$Cdsoil_mgkg)
boxplot(Cacao$Cdsoil_mgkg)
plot(density(Cacao$Cdsoil_mgkg))
plot(log(Cacao$Cdsoil_mgkg))

Revisión del Modelo de Elevación Digital:

DEM <- raster("data_IIC1/DEM_DRONE/Zona_Oeste.tif")
str(DEM)
## Formal class 'RasterLayer' [package "raster"] with 12 slots
##   ..@ file    :Formal class '.RasterFile' [package "raster"] with 13 slots
##   .. .. ..@ name        : chr "D:\\3.0Adicionales\\FIGMM_SMAUML\\ParteII\\SlidesR3\\data_IIC1\\DEM_DRONE\\Zona_Oeste.tif"
##   .. .. ..@ datanotation: chr "INT1U"
##   .. .. ..@ byteorder   : chr "little"
##   .. .. ..@ nodatavalue : num -Inf
##   .. .. ..@ NAchanged   : logi FALSE
##   .. .. ..@ nbands      : int 4
##   .. .. ..@ bandorder   : chr "BIL"
##   .. .. ..@ offset      : int 0
##   .. .. ..@ toptobottom : logi TRUE
##   .. .. ..@ blockrows   : int [1:4] 256 256 256 256
##   .. .. ..@ blockcols   : int [1:4] 256 256 256 256
##   .. .. ..@ driver      : chr "gdal"
##   .. .. ..@ open        : logi FALSE
##   ..@ data    :Formal class '.SingleLayerData' [package "raster"] with 13 slots
##   .. .. ..@ values    : logi(0) 
##   .. .. ..@ offset    : num 0
##   .. .. ..@ gain      : num 1
##   .. .. ..@ inmemory  : logi FALSE
##   .. .. ..@ fromdisk  : logi TRUE
##   .. .. ..@ isfactor  : logi FALSE
##   .. .. ..@ attributes: list()
##   .. .. ..@ haveminmax: logi FALSE
##   .. .. ..@ min       : logi NA
##   .. .. ..@ max       : logi NA
##   .. .. ..@ band      : int 1
##   .. .. ..@ unit      : chr ""
##   .. .. ..@ names     : chr ""
##   ..@ legend  :Formal class '.RasterLegend' [package "raster"] with 5 slots
##   .. .. ..@ type      : chr(0) 
##   .. .. ..@ values    : logi(0) 
##   .. .. ..@ color     : logi(0) 
##   .. .. ..@ names     : logi(0) 
##   .. .. ..@ colortable: logi(0) 
##   ..@ title   : chr(0) 
##   ..@ extent  :Formal class 'Extent' [package "raster"] with 4 slots
##   .. .. ..@ xmin: num 518578
##   .. .. ..@ xmax: num 519198
##   .. .. ..@ ymin: num 9031751
##   .. .. ..@ ymax: num 9032476
##   ..@ rotated : logi FALSE
##   ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
##   .. .. ..@ geotrans: num(0) 
##   .. .. ..@ transfun:function ()  
##   ..@ ncols   : int 24794
##   ..@ nrows   : int 28989
##   ..@ crs     :Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr "+proj=utm +zone=18 +south +datum=WGS84 +units=m +no_defs"
##   .. .. ..$ comment: chr "PROJCRS[\"WGS 84 / UTM zone 18S\",\n    BASEGEOGCRS[\"WGS 84\",\n        DATUM[\"World Geodetic System 1984\",\"| __truncated__
##   ..@ history : list()
##   ..@ z       : list()
# DEM
# slope <- terrain(DEM, opt ="slope")
# aspect <- terrain(DEM, opt="aspect")
# DEM.hill <- hillShade(slope, aspect, angle = 40, direction = 270)
# esta es importante pero demora       DEM <- crop(DEM, extent(518600, 519200, 9031800, 9032400))
#plot(DEM, main = "Drone - Digital Elevation Model (DEM)", col = topo.colors(20),
#     zlim=c(0, 250))
image(DEM, main = "Drone - Digital Elevation Model (DEM)", col = topo.colors(20),
      zlim=c(0, 250))

Revisión del DEM con factores:

Factores_cacao <- st_read(dsn = "data_IIC1/Cacao_Factors/Age_Planta.shp")
## Reading layer `Age_Planta' from data source 
##   `D:\3.0Adicionales\FIGMM_SMAUML\ParteII\SlidesR3\data_IIC1\Cacao_Factors\Age_Planta.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 6 features and 5 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 518541 ymin: 9031841 xmax: 519584 ymax: 9032609
## Projected CRS: WGS 84 / UTM zone 18S
str(Factores_cacao)
## Classes 'sf' and 'data.frame':   6 obs. of  6 variables:
##  $ Id        : num  0 1 2 3 4 5
##  $ Age       : chr  "E1" "E5" "E6" "E4" ...
##  $ Shape_Leng: num  1085 1103 806 1199 1108 ...
##  $ Shape_Area: num  69566 73596 38867 66861 51878 ...
##  $ Uso_Tierra: chr  "Pasto 10" "Descanso 3" "Maiz 10" "Maiz 10" ...
##  $ geometry  :sfc_POLYGON of length 6; first list element: List of 1
##   ..$ : num [1:7, 1:2] 519345 519313 519527 519584 519525 ...
##   ..- attr(*, "class")= chr [1:3] "XY" "POLYGON" "sfg"
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA
##   ..- attr(*, "names")= chr [1:5] "Id" "Age" "Shape_Leng" "Shape_Area" ...
summary(Factores_cacao)
##        Id           Age              Shape_Leng       Shape_Area   
##  Min.   :0.00   Length:6           Min.   : 805.5   Min.   :38867  
##  1st Qu.:1.25   Class :character   1st Qu.:1084.5   1st Qu.:52966  
##  Median :2.50   Mode  :character   Median :1093.8   Median :61545  
##  Mean   :2.50                      Mean   :1064.0   Mean   :59500  
##  3rd Qu.:3.75                      3rd Qu.:1106.5   3rd Qu.:68890  
##  Max.   :5.00                      Max.   :1198.6   Max.   :73596  
##   Uso_Tierra                 geometry
##  Length:6           POLYGON      :6  
##  Class :character   epsg:32718   :0  
##  Mode  :character   +proj=utm ...:0  
##                                      
##                                      
## 
Factores_cacao <- Factores_cacao %>% mutate_if(is.character, as.factor)
summary(Factores_cacao)
##        Id       Age      Shape_Leng       Shape_Area                Uso_Tierra
##  Min.   :0.00   E1:1   Min.   : 805.5   Min.   :38867   Descanso 3       :1   
##  1st Qu.:1.25   E2:1   1st Qu.:1084.5   1st Qu.:52966   Maiz 10          :3   
##  Median :2.50   E3:1   Median :1093.8   Median :61545   Maiz 7 Descanso 3:1   
##  Mean   :2.50   E4:1   Mean   :1064.0   Mean   :59500   Pasto 10         :1   
##  3rd Qu.:3.75   E5:1   3rd Qu.:1106.5   3rd Qu.:68890                         
##  Max.   :5.00   E6:1   Max.   :1198.6   Max.   :73596                         
##           geometry
##  POLYGON      :6  
##  epsg:32718   :0  
##  +proj=utm ...:0  
##                   
##                   
## 
mapview(Factores_cacao)
mapview(DEM, legend =TRUE)+
  mapview(Factores_cacao, alpha.regions = 0.01 , legend.opacity = 0.10, 
         lwd = 2, color = "blue")
## Warning in rasterCheckSize(x, maxpixels = maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ; 
## the supplied Raster* has 718753266 
##  ... decreasing Raster* resolution to 5e+05 pixels
##  to view full resolution set 'maxpixels =  718753266 '
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition

Generando Shapefiles (Opcional):

Factores_cacao2 <- as.data.frame(Factores_cacao)
#Factores_cacao2$Uso_Tierra <- edit(Factores_cacao2$Uso_Tierra)

Factores_cacao2$Var_cocoa <- c("Creole-Aromatic","CCN51","","CCN51","Creole-Aromatic",
                               "CCN51")
Factores_cacao3 <- st_as_sf(Factores_cacao2)

mapview(Factores_cacao3, zcol="Var_cocoa")
Honoria_elementos <- readxl::read_xlsx(path = "data_IIC1/BD_Cd.xlsx")
str(Honoria_elementos)
## tibble [24 × 7] (S3: tbl_df/tbl/data.frame)
##  $ Xeast         : num [1:24] 518642 518644 518690 518701 518706 ...
##  $ Ynorth        : num [1:24] 9032291 9032343 9032334 9032232 9032284 ...
##  $ Cdsoil_mgkg   : num [1:24] 1.75 1.6 1.9 1.7 1.25 1.75 1.5 1.4 1.4 1.35 ...
##  $ Pbsoil_mgkg   : num [1:24] 5.41 5.03 7.24 5.81 5.91 5.62 6.57 5.26 5.01 7.02 ...
##  $ Znsoil_mgkg   : num [1:24] 70.5 99.4 77.1 97.1 74.7 ...
##  $ pHapprox      : num [1:24] 7.06 6.55 6.28 8 6.7 6.4 6.53 5.79 5.41 6.39 ...
##  $ Ceapprox_uS_cm: num [1:24] 39.9 10.3 5 126.7 8 ...
Honoria_elementos <- Honoria_elementos %>% rename(Este= Xeast, Norte=Ynorth, 
                                                  Cd = Cdsoil_mgkg,
                                                  Pb = Pbsoil_mgkg,
                                                  Zn = Znsoil_mgkg,
                                                  pH = pHapprox,
                                                  Ce = Ceapprox_uS_cm)
Honoria_elementos2 <- Honoria_elementos[ ,c("Norte","Este")]
Honoria_elementos2 <- Honoria_elementos2[ ,order(c(names(Honoria_elementos2)))]
sputm  <- SpatialPoints(Honoria_elementos2, proj4string=CRS("+proj=utm +zone=18 +south +datum=WGS84")) 
spgeo  <- spTransform(sputm, CRS("+proj=longlat +datum=WGS84"))
spgeo  <- as.data.frame(spgeo)
colnames(spgeo) <- c("lng","lat")
Honoria_elementos2 <- cbind(Honoria_elementos,spgeo)

Honoria_elementos3 <- st_as_sf(Honoria_elementos2, coords = c("lng", "lat"),
                               remove = FALSE, crs = 4326, agr = "constant")

mapview(Honoria_elementos3, zcol = "Cd") +mapview(Factores_cacao3, zcol="Var_cocoa")+
  mapview(Factores_cacao3, zcol = "Uso_Tierra")+
  mapview(Factores_cacao3, zcol = "Age")
Factores_cacao4 <- Factores_cacao3 %>%  st_transform(crs = 4326 )
  
H <- Honoria_elementos3 %>%  st_intersection(Factores_cacao4)

mapview(H, zcol = "Uso_Tierra", cex = "Cd")+
  mapview(Factores_cacao, alpha.regions = 0.01 , legend.opacity = 0.10, 
          lwd = 2, color = "blue")
mapview(H, zcol = "Age", cex = "Cd") +
  mapview(Factores_cacao, alpha.regions = 0.01 , legend.opacity = 0.10, lwd = 2, color = "blue")
mapview(H, zcol = "Var_cocoa", cex = "Cd") +
  mapview(Factores_cacao, alpha.regions = 0.01 , legend.opacity = 0.10, lwd = 2, color = "blue")

Estandar de Calidad de Cd en Suelos:

Valor límite del Cd en suelos es: \[Cd<=1.4\]

mapview(Honoria_elementos3, zcol = "Cd", at = c(0,1.4,2.5), legend = TRUE,
         col.regions = c("green", "red"))+
  mapview(Factores_cacao, alpha.regions = 0.01 , legend.opacity = 0.10, 
          lwd = 2, color = "blue")

Data Table:

Factores_cacao_n <- Factores_cacao 
Honoria_elementos_n <- Honoria_elementos3 %>% st_transform(crs = st_crs(Factores_cacao))
# sf_use_s2(TRUE) #cambiar a spherical geometry
intersection <- st_intersection(Honoria_elementos_n, Factores_cacao_n)
intersection
Honoria_statical <- intersection %>% st_set_geometry(NULL)

Statical Summary:

min_max_mean_sd <- list(
  min = ~min(.x, na.rm = TRUE),
  max = ~max(.x, na.rm = TRUE),
  mean = ~mean(.x, na.rm = TRUE),
  sd = ~sd(.x, na.rm = TRUE)
)
Honoria_statical %>% 
  summarise(across(c(Cd, Pb, Zn, pH, Ce), min_max_mean_sd))
Honoria_statical %>% group_by(Age) %>% 
  summarise(across(c(Cd, Pb, Zn, pH, Ce), min_max_mean_sd))

Ver vignette("colwise") para más detalles.

Boxplot Analysis:

Cd total:

ggplot(data = Honoria_statical)+
  geom_boxplot(mapping = aes(x = Cd))+
  coord_flip()

Cd by Factor Age and Land Use:

ggplot(data = Honoria_statical)+
  geom_boxplot(mapping = aes(x = Cd))+
  facet_grid(. ~ Age)+
  coord_flip()

ggplot(data = Honoria_statical)+
  geom_boxplot(mapping = aes(x = Cd))+
  facet_grid(. ~ Uso_Tierra)+
  coord_flip()

Cd by Prior Exploration Field:

ggplot(data = Honoria_statical)+
  geom_boxplot(mapping = aes(x = Cd))+
  facet_grid(. ~ Uso_Tierra)+
  coord_flip()

Average Cadmiunm Soil

Corrplot Bivariant

correlacion <- cor(Honoria_statical[ ,c("Cd", "Pb", "Zn", "pH", "Ce")]
                   , method = "pearson")
corrplot(correlacion, method ="circle", diag = TRUE, title= "Correlacion de Pearson Honoria", 
         hclust.method = "median", addrect = 2)
corrplot(correlacion, method ="circle", diag = TRUE, title= "Correlacion de Pearson Honoria")

correlacion2 <- cor(Honoria_statical[ Honoria_statical$Uso_Tierra=="Maiz 10", c("Cd", "Pb", "Zn", "pH", "Ce")], method = "pearson")
# correlacion2 <- cor(Honoria_statical %>% filter(Uso_Tierra=="Maiz 10") %>%
#                     select(Cd, Pb, Zn, pH, Ce), method = "pearson")
corrplot(correlacion2, method ="number")

corrplot(correlacion2, method ="circle", order = "AOE")

corrplot(correlacion2, method = 'shade', order = 'AOE', diag = FALSE) 

corrplot(correlacion2, method = 'square', order = 'AOE', diag = FALSE, type = "lower")

correlacion3 <- cor(Honoria_statical[Honoria_statical$Uso_Tierra=="Pasto 10", c("Cd", "Pb", "Zn", "pH", "Ce")], method = "pearson")
corrplot(correlacion3, method ="circle")

corrplot.mixed(correlacion3, order = 'AOE')

Para ver toda la funcionalidad revisar : An Introduction to corrplot Package.

Opcional

Biplot Multivariant (Principal Component Analysis)

Cadmium Final Geostatistics Map:

library(gstat)
class(Honoria_elementos3)
## [1] "sf"         "data.frame"
summary(Honoria_elementos3)
##       Este            Norte               Cd              Pb       
##  Min.   :518642   Min.   :9031868   Min.   :1.100   Min.   :4.740  
##  1st Qu.:519033   1st Qu.:9031913   1st Qu.:1.350   1st Qu.:5.247  
##  Median :519183   Median :9032036   Median :1.575   Median :5.850  
##  Mean   :519108   Mean   :9032052   Mean   :1.634   Mean   :5.929  
##  3rd Qu.:519221   3rd Qu.:9032134   3rd Qu.:1.863   3rd Qu.:6.525  
##  Max.   :519556   Max.   :9032343   Max.   :2.550   Max.   :7.240  
##        Zn              pH              Ce              lng        
##  Min.   :11.80   Min.   :5.410   Min.   :  4.50   Min.   :-74.83  
##  1st Qu.:44.79   1st Qu.:6.160   1st Qu.:  9.00   1st Qu.:-74.83  
##  Median :58.70   Median :6.375   Median : 13.50   Median :-74.83  
##  Mean   :58.66   Mean   :6.494   Mean   : 26.11   Mean   :-74.83  
##  3rd Qu.:71.89   3rd Qu.:6.588   3rd Qu.: 40.02   3rd Qu.:-74.83  
##  Max.   :99.40   Max.   :8.000   Max.   :126.69   Max.   :-74.82  
##       lat                  geometry 
##  Min.   :-8.758   POINT        :24  
##  1st Qu.:-8.758   epsg:4326    : 0  
##  Median :-8.757   +proj=long...: 0  
##  Mean   :-8.757                     
##  3rd Qu.:-8.756                     
##  Max.   :-8.754

Bubbgle Map:

Honoria_elementos4 <- Honoria_elementos3 %>% st_set_geometry(NULL)
coordinates(Honoria_elementos4) <- ~Este+Norte
class(Honoria_elementos4)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
summary(Honoria_elementos4)
## Object of class SpatialPointsDataFrame
## Coordinates:
##           min     max
## Este   518642  519556
## Norte 9031868 9032343
## Is projected: NA 
## proj4string : [NA]
## Number of points: 24
## Data attributes:
##        Cd              Pb              Zn              pH       
##  Min.   :1.100   Min.   :4.740   Min.   :11.80   Min.   :5.410  
##  1st Qu.:1.350   1st Qu.:5.247   1st Qu.:44.79   1st Qu.:6.160  
##  Median :1.575   Median :5.850   Median :58.70   Median :6.375  
##  Mean   :1.634   Mean   :5.929   Mean   :58.66   Mean   :6.494  
##  3rd Qu.:1.863   3rd Qu.:6.525   3rd Qu.:71.89   3rd Qu.:6.588  
##  Max.   :2.550   Max.   :7.240   Max.   :99.40   Max.   :8.000  
##        Ce              lng              lat        
##  Min.   :  4.50   Min.   :-74.83   Min.   :-8.758  
##  1st Qu.:  9.00   1st Qu.:-74.83   1st Qu.:-8.758  
##  Median : 13.50   Median :-74.83   Median :-8.757  
##  Mean   : 26.11   Mean   :-74.83   Mean   :-8.757  
##  3rd Qu.: 40.02   3rd Qu.:-74.83   3rd Qu.:-8.756  
##  Max.   :126.69   Max.   :-74.82   Max.   :-8.754
bubble(Honoria_elementos4, "Cd", col = c("green","green"), main = "Cadmium Concentrations (ppm)")

Grilla generada por Meuse Data:

data(meuse.grid)
ggplot(data = meuse.grid) + geom_point(aes(x, y))

# Hacemos una grilla regular con `SpatialPolygonsDataFrame`
proj4string(Honoria_elementos4) <- CRS("+proj=utm +zone=18 +south +datum=WGS84") #definimos proyección
Honoria_elementos4  <- spTransform(Honoria_elementos4, CRS("+proj=longlat +datum=WGS84"))

grd <- makegrid(Honoria_elementos4, n = 1000) #generamos grilla
colnames(grd) <- c("x", "y") #asignamos nombres
plot(grd$x,grd$y) #visualizamos


# Next, convert the grid to `SpatialPoints` and subset these points by the polygon.
grd_pts <- SpatialPoints(
  coords      = grd, 
  proj4string = CRS(proj4string(Honoria_elementos4))
)


# Then, visualize your clipped grid which can be used for kriging
ggplot(as.data.frame(coordinates(grd_pts))) +
  geom_point(aes(x, y))

gridded(grd_pts) <- TRUE
grd_pts <- as(grd_pts, "SpatialPixels")
Limite2 <- readOGR(dsn="data_IIC1/Cacao_shp//Limite_Cacao..shp")
## OGR data source with driver: ESRI Shapefile 
## Source: "D:\3.0Adicionales\FIGMM_SMAUML\ParteII\SlidesR3\data_IIC1\Cacao_shp\Limite_Cacao..shp", layer: "Limite_Cacao."
## with 1 features
## It has 1 fields
## Integer64 fields read as strings:  Id
grd <- makegrid(Limite2, n = 1000) #generamos grilla
colnames(grd) <- c("x", "y") #asignamos nombres
plot(grd$x,grd$y) #visualizamos

# Next, convert the grid to `SpatialPoints` and subset these points by the polygon.
grd_pts <- SpatialPoints(
  coords      = grd, 
  proj4string = CRS(proj4string(Limite2))
)

# subset all points in `grd_pts` that fall within `spdf`
grd_pts_in <- grd_pts[Limite2, ]

grd_pts_in2 <- grd_pts_in
proj4string(grd_pts_in2) <- CRS("+proj=utm +zone=18 +south +datum=WGS84")
gridded(grd_pts_in2) <- TRUE
grd_pts_in2 <- as(grd_pts_in2, "SpatialPixels")
cd.idw = idw(Cd~1, Honoria_elementos4, grd_pts_in)
## [inverse distance weighted interpolation]
spplot(cd.idw["var1.pred"], main = " Cadmium inverse distance weighted interpolations")

cd.idw = idw(Cd~1, Honoria_elementos4, grd_pts_in2)
## [inverse distance weighted interpolation]
spplot(cd.idw["var1.pred"], main = " Cadmium inverse distance weighted interpolations")

m <- vgm(.59, "Sph", 874, .04) #setear bien
x <- krige(log(Cd)~1, Honoria_elementos4, grd_pts_in2, model = m)
## [using ordinary kriging]
spplot(x["var1.pred"], main = "ordinary kriging predictions")

spplot(x["var1.var"],  main = "ordinary kriging variance")